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^  Abstract 

We  present  a  collection  of  stability  results  for  finite  difference  approximations  to  the  adrection- 
dif fusion  equation^  =»^ux  -fji^u^.  The  results  are  for  centered  difference  schemes  in  space 
and  include  explicit  and  implicit  schemes  in  time  up  to  fourth  order  and  schemes  that  use 
different  space  and  time  discretisations  for  the  advective  and  diffusive  terms.  The  results  are 
derived  from  a  uniform  framework  based  on  the  Schur-Cohn  theory  of  Simple  von  Neumann 
Polynomials  and  are  necessary  and  sufficient  for  the  stability  of  the  Cauchy  problem.  Some  of 
the  results  are  believed  to  be  new.  — 
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1.  Introduction 

The  linear  advection-diffosion  equation: 

ut  —  anx  +  bu^,  b  >  0,  (1) 

is  often  used  as  a  model  equation  in  computational  physics,  partly  because  it  models  two  of  the 
most  basic  processes  in  a  physical  system,  namely  adveetion  and  diffusion.  In  this  paper,  we  are 
interested  in  the  stability  analysis  of  approximation  schemes  for  solving  this  model  equation.  An 
understanding  of  the  stability  properties  of  a  computational  scheme  is  important  for  both 
theoretical  questions  of  convergence  and  for  practical  questions  of  sensitivity  to  round-off  errors. 

Since  stability  results  for  many  common  schemes  for  approximating  the  wave  equation  ut  ■ 
aux  and  the  heat  equation  ut  ■  buM  are  well-known  [11],  an  often  used  practical  strategy  is  to 
take  the  more  restrictive  of  the  two  stability  constraints  for  the  wave  and  heat  equations  as  the 
stability  condition  for  the  full  advection-diffusion  equation  (1).  However,  the  stability  results  for 
schemes  approximating  the  equation  ut  ™  aux  +  buu  cannot  always  be  inferred  from  those  for 
the  wave  and  heat  equations.  Moreover,  there  is  a  danger  of  arriving  at  a  condition  that  is  more 
restrictive  than  necessary.  For  example,  it  is  well  known  that  Euler’s  method  for  the  wave 
equation  is  unconditionally  unstable,  but  the  scheme  applied  to  the  advection-diffusion  equation 
(Scheme  E2E2  in  Section  4)  is  actually  conditionally  stable.  Worse  yet,  one  can  easily  arrive  at  a 
condition  that  is  not  sufficient.  For  example,  the  stability  condition  of  the  scheme  that  consists 
of  the  Leap-Frog  method  applied  to  the  ux  term  and  Euler’s  method  applied  to  the  uu  term 
(Scheme  LF2E2  in  Section  4)  is  actually  more  restrictive  than  those  of  the  corresponding  methods 
applied  to  the  wave  and  heat  equations  separately. 

The  definition  of  stability  that  we  employ  here  is  a  generalisation  of  the  classical  von  Neumann 
stability  condition  and  is  designed  to  guarantee  that  the  computed  solution  inherits  one 
important  property  of  the  exact  solution:  that  its  norm  remains  bounded.  We  used  a  unified 
approach  for  deriving  the  stability  results  which  is  based  on  the  Schur-Cohn  theory  of  locating 
zeros  of  polynomials  in  terms  of  their  coefficients.  We  apply  this  technique  to  analyse  a 
collection  of  commonly  used  finite  difference  schemes  that  includes  higher  order  approximations 
in  both  space  and  time. 

Stability  analysis  for  difference  approximations  to  time  dependent  partial  differential  equations 
is  often  tricky,  tedious  and  difficult.  In  this  regard,  it  may  be  of  interest  to  point  out  here  that 
we  have  found  an  erroneous  stability  result  for  the  Euler  scheme  E2E2  given  originally  by  Fromm 
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(  [3),  p.365)  and  later  quoted  by  Roache  (  [13],  p.44)  and  another  erroneous  result  given  in  Roaehe 
(  [13],  p.fil)  for  the  LF2E2  scheme.  The  Schur-Cohn  technique  that  we  employ  here,  however,  is 
extremely  powerful  and  general  (especially  for  analysing  schemes  that  span  more  than  two  time 
levels)  and  can  be  used  in  a  systematic  way  to  derive  stability  results  for  schemes  (not  necessarily 
finite  difference  schemes)  that  are  not  analysed  here. 

Some  of  the  results  that  we  shall  present  here  are  well-known  and  can  be  found  in  books  such 
as  Richtmyer  and  Morton  [11],  Roache  [13],  Vichnevetsky  and  Bowles  [17].  However,  we  believe 
that  some  of  the  results  are  new.  In  any  case,  we  hope  that  the  collection  of  stability  results  in 
this  paper  will  prove  to  be  a  useful  reference. 

In  Section  2,  we  review  centered  difference  approximations  for  the  advective  and  diffusive 
terms.  The  general  framework  of  stability  analysis  and  the  Schur-Cohn  theory  will  be  presented 
in  Section  3.  Analysis  and  results  for  a  collection  of  commonly  used  schemes  will  be  given  in 
Section  4. 

2.  Centered  Difference  Approximations 

In  this  section,  we  collect  for  reference  purpose,  some  well-known  results  concerning  centered 
finite  difference  operators  for  approximating  the  terms  ux  and  u^.  Define  the  translation 


operator: 

T(h)u(x)  =*  u(x+h),  h  >  0.  (2) 

We  can  now  define  the  following  difference  operators  in  terms  of  T(h): 

D+(h)-(T(h)-T(0))/h,  (3) 

D(h)-(T(0)-T(-h))/h,  (4) 

D^h)  -  (T(h)  -  T(-h))  /  2h  -  (D++D  )  /  2.  (5) 


(Notation:  When  the  argument  of  a  difference  operator  is  left  out,  it  is  understood  to  be  h.) 

Approximations  for  ux  and  uu  using  centered  differences  are  well-known  and  are  contained  in 
the  following  theorems: 

Theorem  1:  Formally,  the  first  derivative  Dx  as  d/dx  has  the  following  expansion: 

Dx  -  D0  E  (-1  y  <r.  (h2D+D./4)i, 
where  i=0  J 

*j-!(j!)2  22i]/(2j+l)!. 


(8) 

(7) 


Proof:  See  Kreiss  and  Oliger  [7],  Fornberg  [2]  and  Vichnevetsky  and  Bowles  [17]. 
Theorem  2:  Formally,  the  second  derivative  m  iP/dx"  has  the  following 
expansion: 


where 


D«“D+D-  £0  h  (h  D+D-/4),> 


/»j-[(j!)2  22M/[(2j+l)!(j+l)]. 


(8) 

(9) 


Proof:  See  Swarts  [16]. 

We  shall  denote  the  2m-th  order  difference  approximation  for  Dx  and  Du  by  A^m  and  B2m 
respectively. 

Definition  3:  For  m  >  1,  define: 

2r 


B2m  *  D+D.  E# (-1)1  Mj  (h"D+D/4)>, 


For  the  stability  analysis,  we  shall  need  the  Fourier  transforms  of  these  operators.  We  shall 
define  the  Foorier  transform  of  an  operator  A  by 
A  =  (A  eiqx)  /  eiqx, 


(10) 


where 


q  »■»  2jtw. 


By  noting  that 


D0eiqx  -  ((ei#-e  i#)/(2h))  elqx  -  (i  sin*  /  h)  e,qx, 


and 

where 


(h2D+D  /4)eiqx  -  ((ei,-2+ei#)/4)  eiqx  -  -(sin*(f/2))  eiqx, 
9  =■  2  nub, 

we  can  easily  derive  the  following: 

“  “I  (9in#  /  0)  °\  (s'I>2(#/2))i, 

B2m  -  -  q2  [sin2(f/2)/(f/2)2]  S  ^  (sin2(f/2))i. 


(11) 

(12) 

(13) 


(14) 

(15) 


We  note  that  A-  is  always  purely  imaginary  and  Bj  is  always  real. 

The  coefficients  tr.  and  #».  in  (7)  and  (9)  are  tabulated  for  orders  up  to  six  (i.e.  j  »  0,  1,  2)  in 
Table  2*1.  For  computational  purposes,  it  is  often  more  convenient  to  transform  equations  (6) 
and  (8)  into  ateneil  forms  as: 

XjT(jh)  |  /  h, 


(10) 


v  v  v  -- 


3.  Stability  Analysis  -  General  Discussions 

8.1.  Definition  of  Stability 
We  shall  consider  the  following  Cauchy  problem  for  (1): 
ut  —  aax  +  bu^,  0  <  x  <  1, 

u(0,t)  —  u(l,t). 

u(x,0)  =—  f(x)  m  E  f(cj)eiqx. 

The  exact  solution  is  given  by: 

u(x,t)  =*  E  f(o;)es(u')teiqx, 


where 


s(c*/)  =*  iaq  -  bq2. 


Thus,  the  wave  with  frequency  w  travels  with  speed  a  and  decays  at  an  exponential  rate  given  by 
e‘bq2. 

We  discretize  the  spatial  interval  by  a  uniform  grid  with  mesh  size  h  “  l/(2n+l)  and  use  k  to 
denote  the  time  step.  The  most  general  difference  approximation  to  the  system  (18)  is  of  the 
form  : 

f !  v(x,t+k)  -  E  ^  v(x,t-jk),  (21) 

1  j=o  * 

where  the  s  ^.(a,b,h,k),  j  =■  -l,0,l,..,p  are  spatial  difference  operators: 

£  ^(a,bIh,k)T*'(h).  (22) 

J  *=-■, 

The  difference  scheme  defined  above  has  a  stencil  that  spans  p+2  time  levels,  and  on  the  time 
level  t-jk  (j  =  -l,0,l,...p),  it  spans  the  mesh  points  from  x-m^h  to  x+M^h.  See  Figure  3-1. 

We  shall  assume  that  always  exists  and  is  bounded,  so  that  (21)  can  be  solved  for 

v(x,t+k).  This  usually  amounts  to  requiring  the  band  matrix  defined  by  the  linear  operator  4>ml 
to  be  nonsingular.  Any  reasonable  difference  approximation  has  this  property.  Also,  in  practice, 
initial  values  have  to  be  supplied  for  the  time  levels  t  *=  0,k,...pk.  These  values  can  be  supplied 
by  using  one-step  (two-level)  schemes  for  starting,  for  example,  but  for  our  analysis  we  shall 
assume  that  these  values  are  obtained  from  the  exact  solution  u(x,t). 

We  look  for  approximate  solutions  to  (21)  of  the  form 

v(x,mk)—  E  f(w)  Rm(l)  eiqx.  (23) 


Figure  3-1:  Stencil  of  a  General  Difference  Scheme 
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Note  that  if  m  «  0  in  (23),  we  shall  end  op  with  the  correct  initial  function  v(x,0)  =■«  E  f(w)e1<ix. 
It  can  be  shown  (  [8],  Ch.  9)  that  v(x,mk)  as  given  by  (23)  will  satisfy  the  difference  equation  (21) 
if  R(0)  satisfies  the  characteristic  equation: 


where 


i.,Rp+1  -  ^  -  0, 

?.-2  2 


is  the  Fourier  transform  of  .  R(0)  is  usually  called  the  amplification  factor  of  the  difference 
scheme. 

The  characteristic  equation  (24)  is  a  polynomial  equation  of  degree  p+1  in  R  and  so  has  p+1 
roots.  Only  one  of  these,  usually  called  the  principal  root,  corresponds  to  the  approximate 
solution  v(x,mk)  that  we  want.  The  other  p  roots  are  usually  called  spurious  roots.  In  practice, 
any  error  introduced  in  the  computation  will  be  propagated  by  all  the  spurious  roots.  Therefore, 
unless  there  are  some  restrictions  on  the  spurious  roots,  these  propagated  errors  may  become 
unbounded  and  overwhelm  the  approximate  solution  that  we  seek.  Since  the  exact  solution  u(x,t) 
to  the  system  has  a  norm  (or  energy)  that  is  decreasing  with  time  (or  at  least  not  growing  with 
time),  it  seems  reasonable  to  ask  the  same  from  the  approximate  solution  v(x,t).  This  is  what 
Richtmyer  and  Morton  [1 1]  referred  to  as  the  practical  stability  criteria.  It  usually  turns  out  that 
this  condition  will  be  satisfied  if  we  restrict  the  time  step  k  appropriately. 


The  left  hand  side  of  (24)  is  usually  called  the  characteristic  polynomial.  Its  coefficients  are 
functions  of  a,  b,  h,  k  and  #.  In  general,  we  can  express  the  characteristic  polynomial  of  a  p+2 
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time  level  scheme  as: 

H(R)  —  a„  +  atR  +  ...  +  »p+lRp+1-  (26) 

We  shall  call  the  p+1  roots  of  H(R)  Rj,R2,....,R  +1,  with  Rj  being  the  principal  root.  It  is  clear 
from  (23)  that  a  necessary  condition  for  the  computed  solution  v  not  to  be  growing  is  : 

|Rjl  <  1  Vj.  (27) 

This  is  usually  called  the  von  Neumann  Stability  Condition  and  polynomials  with  property  (2) 
are  called  von  Neumann  Polynomials. 

The  von  Neumann  Condition  is  also  sufficient  for  non-growing  solutions  for  all  two  time  level 
(p  =  0)  difference  schemes  with  only  one  dependent  variable  [II].  However,  it  is  not  sufficient  in 
general.  The  insufficiency  mainly  arises  from  the  fact  that  when  p  >  0,  condition  (27)  does  not 
exclude  the  case  of  multiple  roots  on  the  unit  circle.  Therefore  we  have  to  modify  the  von 
Neumann  condition  a  little  bit. 

Definition  4:  We  shall  call  polynomials  H(R)  with  the  following  property: 

IRjIci  vj, 

Schur  Polynomials. 

We  shall  call  polynomials  H(R)  with  the  following  property: 

|Rjl<i  Vj, 

and  (28) 

Rj  distinct  on  |R|  »  I. 

Simple  von  Neumann  Polynomials. 

Definition  5:  We  shall  call  a  scheme  stable  if  its  characteristic  polynomial  is  a  Simple 

von  Neumann  Polynomial  V  $  £  [0,  2n]. 

Note  that  this  definition  of  stability  is  necessary  and  sufficient  for  the  computed  solution  to 
not  have  a  growing  norm.  Since  the  roots  R.  are  functions  of  a,  b,  h,  k  and  9,  the  stability 
condition  will  impose  a  restriction  on  the  range  of  values  that  the  first  four  parameters  can  take. 

The  notion  of  stability  defined  here  is  analogous  to  the  notion  of  rero- stability  in  the  theory  of 
difference  methods  for  the  initial  value  problem  in  ordinary  differential  equation  (  [0],  p.33  and 
[6],  p.412  ).  Condition  (28)  is  the  so-called  root-condition  in  that  theory. 

Notice  that  our  definition  of  stability  (that  of  non-growing  solutions)  is  slightly  different  from 
the  definition  of  stability  used  in  Richtmyer  and  Morton  [11]  and  in  particular  the  discussion 
about  the  effects  of  lower  order  terms  on  the  stability  for  the  heat  equation  on  p.195  of  their 
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book  does  not  apply  in  our  case.  Their  definition  of  stability  allows  growth  in  the  solution  and, 
for  diffusion  problems  like  ut  *■  b  uH,  stability  is  practically  unaffected  by  lower  order  terms  like 
a  ux. 

3.2.  The  Schur-Cohn  Theory 

There  is  a  whole  theory,  originating  from  Schur  [14,  15]  that  deals  with  tl 
Neumann  Polynomials.  This  theory,  an  excellent  exposition  of  which  can  l 
J.  J.  H.  Miller  [10],  enables  one  to  determine  conditions  on  the  coefficient; 
polynomial  for  it  to  be  Simple  von  Neumann.  We  shall  present  the  main 
here  and  shall  refer  the  reader  to  the  original  papers  for  more  details. 

Given  a  polynomial 

9 

Mi)  =  aA  +  a.z  +  ...  +  a.z*'  a  E  a.zJ, 
w  1  v  j=o  J 

of  degree  v  (a^  ^  0)  and  having  no  sera  at  the  origin  (aQ  ^  0),  (any  given  polynomial  can  be 
reduced  to  this  case  without  losing  information  about  the  location  of  its  zeros),  one  can  associate 
with  <j>  another  polynomial  >  satisfying  the  same  conditions,  and  defined  by 

**(»)—  £  Vjli’ 

j“0  * 

where  a  denotes  the  complex  conjugate  of  a.  The  reduced  polynomial  <j>1  is  defined  by 

^,(»)  —  (  ^*(0)  fa)  *  *(0)  **(*)  )  /  *•  (29) 

The  main  results  that  we  need  are  contained  in  the  following  two  theorems: 

Theorem  6:  ^  is  a  Schur  Polynomial  iff  |^*(0)|  >  |^(0)|  and  ^  is  a  Schur  Polynomial. 
Theorem  7:  ^  is  a  Simple  von  Neumann  Polynomial  iff  either  |^*(0)|  >  |^(0)|  and 
is  a  Simple  von  Neumann  Polynomial  or  as  0  and  is  a  Schur  Polynomial  (^’ 
denotes  the  derivative  of  ^  with  respect  to  its  dependent  variable). 

By  repeated  applications  of  the  above  two  theorems,  it  is  possible  to  reduce  the  question  of 
whether  a  n-th  degree  polynomial  is  a  Simple  von  Neumann  Polynomial  to  that  for  a  first  degree 
polynomial,  which  can  be  solved  solved  more  easily  by  analytical  means.  These  results  turn  out 
to  be  very  useful  for  determining  stability  limits  of  difference  schemes,  as  compared  to  first 
finding  the  roots  of  the  characteristic  polynomial  explicitly  and  then  determining  their  absolute 
values.  Furthermore,  this  last  approach  may  not  even  be  applicable  for  polynomials  of  higher 
degrees. 


lass  of  Simple  von 
>und  in  a  paper  by 
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4.  Stability  Analysis  -  Specific  Schemes 


4.1.  Some  Commonly  Used  Schemes 

In  this  section,  we  shall  present  the  stability  analysis  and  results  for  some  commonly  used 
difference  schemes  for  solving  the  advection-diffusion  equation.  We  shall  adopt  the  following 
convention  for  naming  the  schemes. 

Definition  8:  The  name  for  a  scheme  shall  consist  of  four  fields: 

Scheme  Name:  A  B  C  D 

where  AB  is  used  to  denote  how  the  scheme  treats  the  advective  term  aux,  and  CD  is 
used  to  denote  how  the  scheme  treats  the  diffusive  term  bu^.  A  and  C  are  used  to 
denote  the  time  discretization  method  used  for  the  aux  term  and  the  bu^  term 
respectively.  B  and  D  are  used  to  denote  the  order  of  the  centered  differencing  used  for 
the  aux  and  buxx  terms  respectively. 

The  following  abbreviations  will  be  used  for  the  time  discretizations: 

E  -  Forward  Euler  (first  order,  two  levels,  explicit) 

BE  -  Backward  Euler  (first  order,  two  levels,  implicit) 

CN  •  Crank-Nicolson  (second  order,  two  levels,  implicit) 

LF  -  Leap-Frog  (second  order,  three  levels,  explicit) 

DF  -  DuFort-Frankel  (second  order,  three  levels,  explicit) 

BD  -  Backward  Differencing  (second  order,  three  levels,  implicit) 

P4  -  Pade  (fourth  order,  two  levels,  implicit) 

For  example,  the  following  scheme: 

(vm+1  -  vm)/k  =  aD0vm  +  bD+D  v™ 

will  be  denoted  by  E2E2  because  Euler’s  method  is  used  to  discretize  in  time  and  the  spatial 
approximations  are  second  order. 

We  shall  analyse  the  following  classes  of  schemes:  EnEj,  BEnBEj,  CNnCNj,  P4nP4j,  BDnBDj, 
LFnCNj,  LFnEj  and  LFnDFj,  where  n  and  j  are  even  nonzero  integers.  This  set  of  schemes  is  by 
no  means  exhaustive  but  is  intended  to  include  most  of  the  commonly  used  schemes.  It  includes 
schemes  that  are  first  order,  second  order  and  fourth  order  in  time;  schemes  that  use  the  same 
order  of  spatial  approximation  for  both  the  aux  and  bu^  terms  and  those  that  use  different 
orders  for  the  two  terms;  schemes  that  use  the  same  temporal  scheme  for  both  terms  and  those 
use  different  temporal  schemes  for  them;  explicit  schemes  and  implicit  schemes;  and  finally  two- 
level  and  three-level  schemes. 
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We  collect  in  Table  4-1  the  exact  definitions  of  these  schemes  and  their  stability  conditions. 
For  reference  purpose,  we  have  also  indicated  the  order  of  the  truncation  errors  for  each  scheme. 
For  more  details  on  the  error  analysis  and  the  stencils  for  these  schemes,  the  reader  is  referred  to 
[!]• 


4.2.  Stability  Analysis 

Next,  we  shall  present  the  stability  analysis  for  the  numerical  schemes  presented  in  Table  4-1. 
We  shall  apply  the  basic  Schur-Cohn  theory  presented  in  Section  3.2  to  the  characteristic 
polynomial  of  each  of  the  schemes.  Only  the  three-level  schemes  make  non-trivial  use  of  this 
theory  and  we  shall  present  only  their  analysis  in  details. 


We  shall  use  the  following  definitions  in  this  section: 

Definition  9:  Define: 

P  =  4kb/h2  , 

a  =  ak/h  , 

7  =*  -  bk  B.  , 

*  -  (An/i). 

We  note  that  all  the  above  quantities  are  real  and  f)  and  7  are  non-negative.  The  indices  n  and  j 
should  be  clear  from  the  context. 


We  shall  also  need  the  following  definitions: 

Definition  10:  Define: 

M  =*  min(n,  j), 

./*•! 

*n~ 

v=0 

i/«-« 

'i-  E,'-' 

where  the  nu's  and  the  trjs  are  defined  in  Section  2.  Specifically,  k2 
11/5;  p2  =»  1,  p4  =*  4/3  and  =  68/45. 


1,  k4  =  5/3  and  ic# 


1)  EnEj 


The  amplification  factor  is  given  by 
R  =■  1  +  ii  -  7. 

We  have  to  find  conditions  on  h  and  k  so  that  |R|  <  1  for  6  6  [*w,  n*].  This  leads  to  the 
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Table  4-1:  Summary  of  Stability  Results  for  Schemes  for  ut  aux  +  bu 


Notation:  h:  space  step,  k:  time  step. 

n,  j  :  positive  even  integers,  M  =  »in[n  ,  j]. 

Ar  :  n-th  order  centered  difference  operator  for  u  . 
Bj  :  j-th  order  centered  difference  operator  for  uM. 
Order(p,  q)  :  Truncation  error  =  Q(hp)  ♦  OCk^). 


Scheme 


Definition 


I  Stabi I ity  Condition 
I  (order  in  a,  order  in  t) 


E2E2 


(v**1  -  vB)/k  =  (aA2  ♦  bB2)v“ 


I  k  <  min[2b/a2.  h2/2b] 
I  order(  2.  2) 


BEnBEj 


(v'*1  -  vB)/k  =  (aA.  ♦  bB.Jv^1 
n  J 


I  Unconditionally  Stable 
I  order(  N,  1) 


I 

I 

I 

-+ 


CNnCNj 


(v-*1  -  vB)/k  = 

(aA.  ♦  bB.) (vB+1  ♦  vB)/2 


I  Unconditionally  Stable 
I  order(  N,  2) 

-+ - 

I  Unconditionally  Stable 

I 

I  order(  N,  4) 


P4nP4j 


(I  -  G/2  ♦  G2/12)i~l  = 
(I  ♦  G/2  ♦  G2/12) v" 
where  G  =  k(aA.  ♦  bB;) 

ra  j 


BDnBDj 


(3/2)(vB+1-vB)/k  -  (l/aXv^-v^J/k 


=  (aA  ♦  bB-) v 
■  J 


■♦1 


I  Unconditionally  Stable 
I  order(  H,  2) 


LFnCNj 


(v«-l  _  »,-1)/(2k)  =  aA#vB 
♦  bB  (v**1  ♦  v,_1)/2 


I  Same  as  that  for  LFn: 

I  n  =  2  :  k  <  h/|a| 

I  n  =  4  :  k  <  0.7287  h/|a| 

I  n  =  6  :  k  <  0.6305  h/|a| 

I  order(  M.  2) 


LF2E2 


(v,M  -  v#-l)/(2k)  =  aA2vB  ♦  bB^*-1 


(ak/h)2  ♦  (4bk/h2)  <  1 
order(  2,  2) 


LFnDFj 


(vB*1-v"”1)/(2k)  =  aA  v*  ♦  bBjm 
-  Orb/h2)(vB*1  -  2vB  ♦  v"*1) 

«?2  =  1 
*4  =  4/3 
i/6  =  68/45 


I  n=2.j=2:  k  <  h/lal  I 
I  n=4, j=2:  k  <  0.5311  h/|a|  I 
I  n=2, j=4:  k  <  0.9685  h/|a|  I 
I  n=4, j=4:  k  <  0.5458  h/|a|  I 
I  Error  =  0(h",k2,ifj(k/h)2)  I 


condition: 


(1  -  if  +  *  <  1. 

It  follows  that  two  ntttatary  conditions  are: 

7  <  2  and  |j|  <  1, 

which  reduces  to 

fo.Jt  <  1  and  |a|«0  <  1. 

Sufficient  conditions,  however,  are  more  difficult  to  derive  analytically,  especially  for  larger 
values  of  n  and  j.  For  the  simplest  case  of  n  ■■  j  —  2,  it  can  be  shown  that  the  necessary  and 
sufficient  conditions  are: 

0  <  2  and  a2  <  0/2  , 

which  reduces  to  (30) 

k  <  min(2b/a2,  h2/2b). 

The  analytical  solution  of  this  problem  is  not  difficult  but  a  bit  tedious  and  can  be  found  in  [l]. 
For  a  geometric  proof,  see  [12].  Results  for  the  general  case  are  not  known. 

Remarks:  The  stability  of  the  method  E2E2  was  studied  by  Roac he  [13]  and  a  twin 
dimensional  version  by  Fromm  [3].  Instead  of  condition  (30),  they  found  lower  bounds  on  the 
spatial  step  site  h  independent  of  the  temporal  step  sise  k,  the  so-called  ceil  Reynolds  Number 
limitation,  which  is  more  restrictive.  Our  results  show  that  h  can  be  as  small  as  we  wish.  As 
long  as  k  is  small  enough,  the  scheme  is  stable.  See  also  Hirt  [5]. 

2)  BEnBEj 

The  amplification  factor  is: 

R=»l/(l-W+7). 

IV  ^ 

It  follows  from  the  definitions  of  A„  and  B.  in  Section  2  that 

»  l 

|R|2  -  1  /  (1  +  ^  <  1. 

Hence  this  scheme  is  unconditionally  stable. 

3)  CNnCNj 

The  amplification  factor  is: 

R  -  (  1  +  16/ 2  -  7/2  )  /  (  1  -  10/2  +  7/2  ). 

It  follows  that 

|R|2  -  «l-l/2)2  +  (I/2I2)  /  ((I+1/2)2  +  d/2)2)  <  1. 

Hence  this  scheme  is  unconditionally  stable. 


4)  P4nP4j 

The  amplification  factor  is: 

R  -  (  1  +  G/2  +  G2/12  )  /  (  1  -  G/2  +  G2/12  ). 

Let  G  -■  A  +  i  1  where  Z  and  I  are  real.  Then  it  follows  that 

|R|2  -  ((1-7/2 +z?  +  M+Zf)  /  ((I+7/2+*)2  +  (6/2+Tf)  <  l. 

Hence  this  scheme  is  unconditionally  stable. 

5)  BDnBDj 

In  the  notation  developed  in  Section  3.2,  the  characteristic  polynomial  is  given  by: 

^(z)  —  (3/2  +  7  -  i£)z2  -  2z  +  1/2  , 
and 

) «“  l/2z2  -  2z  +  (3/2  +  7  +  i  $) . 

We  shall  use  Theorem  7  to  show  that  4(z)  is  a  simple  von  Neumann  polynomial.  It  would  then 
follow  that  the  scheme  is  stable. 

The  condition  |^*(0)|  >  |d(0)|  is  certainly  always  satisfied.  We  next  compute  u 
*,(*)  -  [(3/2  +  7)2  +  62  -  1/4]*  -  2(1  +  7  +  i *)  . 

^j(z)  is  simple  von  Neumann  iff 

|2(1  +7  +  i6) |2  <  [(3/2  +  i?  +  J2  -  1/4]2  . 

This  can  easily  be  shown  to  be  true  for  any  real  7  and  6.  Thus  ^(z)  is  simple  von  Neumann  and 
the  scheme  is  unconditionally  stable. 

6)  LFnCNj 

The  characteristic  polynomial  is  given  by: 

4(*)  —  (1  +  7)*2  -  2i 61  -  (1  -  7). 

We  thus  get: 

^*(z)  “  (7  -  l)z2  +  2i6i  +  (1  +  7). 

The  condition  \4  (0)|  >  |^(0)|  reduces  to  |1  +  7  >  |1  -  7|,  which  is  always  true  because  7  is 
positive. 

Next  we  compute  ^t(z)  as 

^,(z)  —  4tz  -  4i7#  . 


\6\<l  .  (31) 

Note  that  condition  (31)  is  the  same  as  the  stability  criterion  for  the  Leap-Frog  scheme  applied  to 
nt  **  aux,  with  n-th  order  centered  differencing  in  space.  The  criteria  are  all  of  the  form  |ak/h| 
<  cn>  where  the  constants  cB  can  be  found  in  Fomberg  [2],  for  example.  In  particular,  for  c2  — 
1,  c4  —  0.7287  and  c#  -«  0.6305. 

7)  LFnEj 

The  characteristic  polynomial  is  given  by: 

*(*)  -  s2  -  2\6i  -  (1  -  2-r). 

We  thus  get: 

*  (27  *  l)*2  +  2i$s  +  1. 

Now  the  condition  |^*(0)|  >  |^(0)|  reduces  to  1  >  |1  -  2-y|  which  will  be  satisfied  iff 

7  <  1  , 

which  reduces  to 

<  1. 

This  is  the  same  as  the  stability  condition  for  Euler’s  method  applied  to  the  ujx  term  with  a  time 
step  of  2k  and  is  clearly  a  necessary  condition  for  stability. 

We  next  compute  ^(z)  as 

^(*)  —  *(1  *  (1  -  27^)  *  4i*7. 

^,(z)  will  be  simple  von  Neumann  iff 

|4$7|  <  1  -  (I-27)2, 

which  reduces  to 

1*1  <  1  -  7- 

This  is  a  relationship  involving  (3,  a  and  #  and  we  want  to  derive  conditions  involving  fi  and  0 
for  it  to  hold  for  all  values  of  #.  In  the  case  of  n  —  2  and  j  ■»  2,  this  has  been  worked  out  in  [1] 
and  the  necessary  and  sufficient  condition  is: 

|a|2  +  0  <  1.  (32) 

Results  for  more  general  values  of  n  and  j  are  not  known. 

Remark:  Roache  [13]  considered  this  scheme  for  n  —  j  »■  2,  but  he  claimed  that  the  stability 
analysis  of  the  advection  and  diffusion  terms  may  be  analysed  separately,  and  thus  obtaining  the 
conditions  |a|  <  1  and  fi  <  1.  (In  his  book,  he  had  the  equivalent  of  $  <  2,  which  I  believe  is 
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either  a  typo  or  a  mere  oversight.)  These  two  separate  conditions  are  much  less  restrictive  than 
(32)  and  we  believe  that  Roache’s  results  were  erroneous.  Riga]  (12)  derived  the  correct  stability 
limit  using  a  different  approach  but  he  did  not  show  that  it  is  both  necessary  and  sufficient. 

8)  LFnDFj 

Generalised  Du  Fort-Frankel  methods  have  been  studied  by  Gottlieb  and  Gustafsson  [4]. 
When  adapted  to  our  case  for  the  equation  ut  -  auz  +  b«  ,  the  scheme  LFnDFj  becomes 
(vm+l  .  vm‘*)/2lc  -  (aAn  +  bBjJv"1  -  (q.b/h2Xvm+l  -  2v“  +  v”*1) 
where  »/■  is  a  positive  constant  chosen  to  make  the  scheme  unconditionally  stable  for  the  case  a 
—  0  (the  heat  equation).  The  conditions  is  exactly: 

>  Pi  •  (33) 

The  truncation  error  of  the  scheme  is  0(hM,  k2,  ^.(k/h)2).  Thus  the  larger  the  value  of  n-}  is,  the 
larger  is  the  truncation  error.  Therefore,  in  what  follows,  we  shall  assume  that  sy.  takes  on  the 
value  of  the  lower  bounds  given  in  (33). 

The  characteristic  polynomial  is: 

p{i)  —  i\\  +  tf.^/2)  -  tin-0  +  2itf  -  2-y)  -  (1  -  q.0/2)  . 

It  follows  that: 

*(s)  —  s^-l  +  q.0/2)  -  s(q.0  -  2\6  -  27)  +  (1  -  q^/2)  . 

The  condition  |^*(0)|  >  |4{0)|  reduces  to 
|1  +  q.0/2|  >  |1  -  f.J/2|  , 

which  is  always  satisfied  because  tfj  and  0  are  both  positive. 

We  compute  ^{(z)  as: 

*,(*)  —  2q.0i  -  2(q.0  -  27  +  nfiS)  . 

Hence  the  stability  criterion  is  given  by 
|1  +  2t/{0n})  +  i$|  <  1  , 

which  can  be  written  as: 

o2  <  (1  -  (1  +  Bj(h2/2,j))2)/|Anh|2  .  (34) 

It  can  easily  be  verified  that,  for  a  given  n  and  j,  the  right  hand  side  of  (34)  is  a  function  of  9 
only,  and  thus  its  minimum  can  be  found,  at  least  numerically,  to  yield  an  upper  bound  for  a2  as 
the  stability  criterion  for  the  scheme  LFnDFj.  Moreover,  it  can  also  be  seen  from  the  form  of 
(34)  that  the  limitation  on  k  is  more  restrictive  than  the  corresponding  limitation  for  the  Leap- 
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Prog  method  applied  to  ut  ■»  aux.  The  upper  bounds  hare  been  competed  numerically  and  their 
values  are  given  in  Table  4-2. 

Table  4-2:  Stability  Constants  for  Scheme  LFnDFj 
Stability  Condition  :  a2  <  c  •  ,  where  c.  .  is  given  below: 

j 


♦- 

1 

2 

- + — 

1 

4 

1 

6 

1 

2  1 

1.0 

1 

0.9686 

1 

0.9242 

1 

4  1 

0.5310 

1 

0.5458 

1 

0.5391 

1 

6  1 

0.3976 

1 

— 4- 

0.4231 

1 

- 

0.4281 

1 
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